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A numerical integration procedure to calculate the second virial coefficient of simple polyatomic molecules is 
proposed. The intermolecular pair potential is assumed to be a sum of a spherically symmetric contribution and 
an unsymmetric, angular dependent, contribution. The method is based on evaluating the possible different 
numerical values for this latter term. Quadrupolar and dipolar molecules are considered. Calculations for the 
virial coefficient for quadrupolar molecules are judged to be correct to within one part in 2500 or better, and to 
within one part in 300 or better for polar molecules. Results from the method are compared to corresponding 
results from the well-known Pople expansion procedure. It is shown that care must be taken to ensure this latter 
technique yields a convergent answer. 
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1 . Introduction 

In several papers [1, 2, 3], 1 we have determined the 
thermophysical properties of simple polyatomic gases via 
statistical mechanics. In particular, a modification of the 
Pople expansion procedure [4] was used to evaluate the 
second virial coefficient. This technique is standard but 
there are some drawbacks when it is applied to molecules 
with relatively large electrostatic moments. The object of 
this paper is to introduce an alternative procedure. 

The second virial coefficient, B, is defined by the expan- 
sion 



where N is Avogadro's number, /3 = 1/kT (k is Boltzmann's 
constant), <f) is the pair potential, and R! and R 2 denote the 
position and angular configuration of molecule 1 and 2, 
respectively. 

In this work the intermolecular interaction of two polya- 
tomic molecules is considered in the simple, but standard 
coordinate system of figure 1: the molecules are treated as 



-^- = 1 + Bp + 
pRT H 



(1) 



R the 



(where p is the pressure, T the temperature, K the gas 
constant and p the density (mol/L)) and is given by an 
expression of the form 



B ~ N jj {exp [-/8<MR 2 , R 2 )] - 1} dRi <*R 2 





* 



(2) 
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Figure 1. The coordinate system for two polyatomic molecules. 

rigid bodies interacting along the line of their centers, dis- 
tance r apart, with a relative orientation given by the angles 
#!, 62 and 4>. The potential can be written as a sum: 
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4>*(r*, U 6 2 , *) = 0*(r*) + 0*(r*, $ u O t , <t>) (3) 

with $*(r*) the spherically symmetric, angular independent, 
part and <f>%g the unsymmetric angular dependent part. In eq 
(3), r* = r/cr and <£>* = </>/€, where (X and € are a length 
and energy parameter of the spherical potential, respectively 
(to be defined below). 

For the coordinate system of figure 1 and the potential of 
eq (3), eq (2) becomes 



B* = B/b = 

- 3/(8tt) sin 0! dd 1 i 
Jo Jo 



do dSi 



f% 



d<& 



X J {exp [- <£*(r*, l9 2 , 4>)/r*] - 1} r* 2 dr* (4) 

where b = (2/3) ttNct 3 and 7* = T/(e/k). 

In previous studies [1-3], we have taken <\>f to be the m- 
6-8 function of Klein and Hanley [5] 



*V) = ^ (' 

m — 6 \ 



m 



y(ro 



8) /o-\ 6 



m — 6 



(r)* 



r(')V (5) 



where d = r min /cr, with cr and r min , and also the energy €, 
defined by the conditions </>(o") = and </>(r min ) = — e. The 
parameters m and y represent the "hardness" of the repulsive 
term, and the "strength" of the inverse eight attractive term, 
respectively. 

If we consider non-polarizable axially symmetric quadru- 
polar molecules, <f>%s is given by 



*&(r* 0i, e t , *) = \-^AQi, o* *) 



(6) 



vhere 



/(0 l9 2 , <D) = 

1 — 5cos 2 0x — 5 cos 2 2 — 15 cos 2 x cos 2 2 

+ 2(sin0x sin0 2 cos<l> - 4 cos0x cos0 2 ) (7) 

and 0* 2 is the reduced quadrupole moment (squared); 0* 2 
= 6 2 /(€C7 5 ). 

To be consistent with our earlier calculations, and those 
of other authors, we will first consider the potential given by 
eqs (3), (5), and (6) but with y = and m = 12 for which 
eq (5) becomes the Lennard- Jones, 12-6, potential. 

2. Procedure 

The evaluation of the second virial is carried out in four 
steps as follows. First, the quantity /(0i, 2 , <&) of eq (7) is 
considered. Specifically, the molecules are set in a selected 
initial configuration, e.g., 6 X = 2 = $ = 0, then each of 
the angles is varied incrementally (e.g., by 7r/4) and the 
numerical value of / noted after each incremental step. 
Because of symmetry, many of the values are identical and 



it is only necessary to record, say, / different values^ (i — 
1, /), and their number of occurrences or weight, wtj. In 
general fi and its weight, wt^, are associated with the initial 
configuration C x and the step 7r/n where x = or 7r/2n 
with /i = 2, 4, 8. . . . That the number / is relatively small 
compared to the total number of configurations is shown in 
table 1. 

Table 1. Values of I for various incremental steps and initial 
configurations of two molecules in the coordinate system of figure 1 . 

(a) Quadrupolar interactions, (b) dipolar interactions. 



Incremental 
step it jn 



Number of different (/) values for 
the angle function at two selected 
initial configurations: Cq and C^n 



Number of 

possible values 

for the angle 





^•o ^7r/2n 


luncuon 


(a) Quadrupolar 


7T/4 

tt/8 
tt/16 


11 12 

74 80 

548 576 


128 
1024 
8192 


(b) Dipolar 


tt/4 
tt/8 
tt/16 


13 24 

123 160 

1011 1152 


128 
1024 
8192 



Second, an intergral Ii is defined by 

h= \ exp [-(f)?/T* - ljr* 2 dr* i = 1, / (8) 

to correspond to the r* integral in eq (4). Here <f>i* is the 
potential of eq (3) for a given quadrupole moment associated 
with a particular/^ (0 l9 2 , <I>). 

The integral is evaluated numerically at a selected 7 1 * by 
the one-dimensional Chebyshev integration technique used 
in our previous work [1, 2], 

Third, the value of the reduced second virial coefficient 
for an initial configuration, C x , and an incremental step 
follows from the expression 



B*{T*, e* 2 ) = -3 l ~ 

2 a 



(9) 



flj = (sin 0j sin 2 wt)j 



= 1, /. 



(10) 



Fourth, and last, since 5* of eq (9) is associated with C x 
and the incremental step, a final estimate of the virial 
coefficient is obtained from the results for different steps at 
a particular C x , and then from a repeat of the procedure for 
a different C x . 

2.1 . Variation of I. With the Angle Function 

The calculation outlined above is potentially time-consum- 
ing, despite the simplification due to symmetry, because the 
integral /$ apparently needs to be evaluated / times for each 
step at each configuration. We have, however, observed a 
feature which allows the calculation to be manageable: 
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namely, that I t appears to vary smoothly with the angle 
funetion at a given 7 1 * and 0* 2 . This is illustrated hy figure 
2 in whieh 3/ t is plotted versus f t (0 u 2 , <E>) for the 12-6 
potential with 0* 2 = 0.7 at three representative tempera- 
tures. Hence for a given spherical potential with a given 
G* 2 at 7 1 *, one only has to evaluate a few I t directly: the 
remainder can be estimated to a high precision by interpola- 
tion. 



QUADRUPOLAR 



3Ij -4- 




-8- 



-12 - 



fj(0, 2 <I>) 

FIGURE 2. Plot of 3l { eq (8), versus U^\^),from eq (7), at T* = 0.7, 
2.0, and 20.0 for a 12-6 quadrupolar molecule with 6* 2 = 0.7. 



3. Results 
3.1 . Quadrupolar Molecules 

Sample detailed results are presented in table 2 which 
displays values of Z?* for several incremental steps and initial 
configurations at four reduced temperatures; 0.7, 1.0, 2.0, 
and 20.0. The spherical potential was taken to be the 12-6 
and the reduced quadrupole moment (squared) was set at 
0.2 — which corresponds to the moment (squared) for a 
molecule similar to nitrogen. Our estimated convergent 
values are also shown and are judged to be precise to within 
one part in 2500 at 7* = 0.7 and to one part in 10,000 at 
T* = 20.0. 

A check on the absolute values of the virials listed in 
table 2 was made by comparing them to corresponding 
values from the expansion technique [4]: the expansion 
procedure can be expected to yield essentially an exact 



result for/?* provided the quadrupole moment is small. For 
the potential of eq (3) with equations (5) and (6) one can 
derive [2, 6]: 



639 e* 8 



B*(T*) ~- 


= B* 


21 e* 4 216 e* 6 

4608 O* 10 
4235 7* 5 28 



245 T* 4 



(11) 



TABLE 2. The second virial coefficient of a quadrupolar molecule, d* 2 = 
0.2, as a function of initial configuration and incremental step (see text). 

Comparisons with equivalent values from the expansion of equation (11) are given. 





Incre- 
mental 
step 
(ir/n) 


Initial configuration 


Estimated* 

convergent 
value 


Equation 

(11) 


7 1 * 


Co 


W/2n 


0.7 

1.0 

2.0 

20.0 


tt/4 
tt/8 
tt/16 

tt/4 
tt/8 
tt/16 

tt/4 
tt/8 
tt/16 

tt/4 
tt/8 
tt/16 


-4.80711 
-4.84904 
-4.85425 

-2.57073 
-2.58929 
-2.59140 

-.63202 
-.63775 
-.63828 

.52567 
.52522 
.525193 


-4.86570 
-4.85841 
-4.85637 

-2.59522 
-2.59293 
-2.59220 

-.63862 
-.63857 
-.63845 

.52534 
.52519 
.525189 


-4.8557 
-2.5919 
-0.63840 
.525190 


-4.8556 
-2.5919 
-0.63834 
.525185 



* Estimated by plotting B* versus tt/h for both C and C^n- The 
intersection point of the two curves was taken as the convergent value. 



where Bf is the angular-independent, classical second virial 
at 7 1 *. J x are dimensionless integrals given by 



./, 



f 

Jo 



exp [-</>* (r*)/7*]r* ( 



di*. 



(12) 



Results from eq (11), for the 12-6 potential and 0* 2 = 
0.2, are given in the last column of table 2. Since the 
agreement between corresponding second virial coefficients 
from our procedure and from equation (11) is very close, we 
conclude our procedure is satisfactory. 



3.2. Polar Molecules 

The integration procedure is not, of course, restricted to 
quadrupolar molecules. It was straightforward to extend it to 
polar molecules interacting with a potential of the form 



0*G*, e u e 2 , <t>) = 0? (i*) 



v> 



*2 



g(0 u e 2 ,<!>) (13) 



yhere 



g(0 u 2 , *) = 2 cos0! cos0 2 - sin0! sin0 2 cos<£. (14) 
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/x* is the reduced dipole moment given by [jl* 2 — ^/(ecr 3 ), 
where fi is the dipole moment. $*(t*) is the spherical 
potential of eq (5), as before. 

Calculations followed those for quadrupolar molecules 
using the two significant time-saving arguments, i.e., (a) the 
symmetry of the coordinate systems which reduces substan- 
tially the number of independent values of g{Q\, 2 , <P) 
which have to be considered (see table 1), and (b) that a 
particular value of g, gi, is a smooth function of the r- 
integral, l\ [ eq (8)], see figure 3. Typical results are shown 
in table 3 for three temperatures. The results are those for a 
12-6 potential for ^* 2 = 0.8485 and 2.8284, respectively. 
We have listed the parameter r*, used in Hirschfelder, 
Curtiss, and Bird [7] and elsewhere [t* = 8 _1/2 /x* 2 ]. For 
these dipole moments (squared), T* = 0.3 and 1.0. 



3Ii 




-2.0 



Figure 3. Plot of3I t eq (8), versus g/fl^O), from eq (14), at T* = 
0.7 and T* = 1.0 (solid line), T* =0.7 and r* = 0.3, (dashed line) and 
T* = 2.0 and r* = 0.3 (dotted line). 

TABLE 3. Comparison between the sound virial coefficient for a polar 

molecule calculated by numerical integration, and by the expansion of 

equation (15) 



J 1 * 


B*, integration 


5*, expansion 


t* = 0.3 


0.7 

2.0 

20.0 


-6.001 
-0.7130 
.5241 


-6.0182 
-0.7168 

.52441 


r* = 1.0 


0.7 

2.0 

20.0 


-41.61 
-1.706 
.5148 


-39.89 
-1.7264 
.51468 



Again our results can be compared to those from the 
expansion technique [6]: 



B*(T*) 

=^- i9 4[w^] w 4- 



>2(j - 3)-| 



J*2j 



yhere 



Kj= ^{2k)!/{k!f 

fc=0 



(15) 



(16) 



and J 6J = J x of eq (12). Values of B* from eq (15) fory = 4 
are given in table 3. 



The results from our integration procedure are apparently 
low. But this is most probably due to the appreciable 
contribution of the term (n* 2 lr 2 )g(6 u 6 2 , 4>) to the integral 
eq (8) at large values of r. Our Chebyshev integration 
method can fail to yield a fully convergent answer in such 
cases. Clearly, one could compensate for this difficulty. 
However, the substantial discrepancy at 7 1 * = 0.7 for r* = 
1.0 is undoubtedly a failure of the expansion (15) to give a 
convergent fory = 4. 

3.3. Mixed Moments 

The method can be applied to molecules with more than 
one electrostatic moment and, in fact, has been used for 
nitrous oxide: a molecule which has a small dipole and a 
large quadrupole [3]. We considered the potential 

0V, 0, , e 2 , *) = <# + <tf, D + 4>% Q + <I>%d (17) 

where (f)f is eq (5); 4>$d is the angle dependent term of eq 
(13); <fi%Q is given by eq (6); while Ofaj), the contribution of 
the quadrupole-dipole interaction, is given by 



<f>%D 



3 

= I7i LC *Q*[( COS ^2 — COS^i)(3 COS#! COS0 2 

2r* 

— 2 sind 1 sin9 2 cos<$>) — cosd x + cos# 2 ]- 



(18) 



We will not report here results for the second virial 
coefficient of nitrous oxide; suffice to say that agreement 
between our method and that from the corresponding expan- 
sion [6] was within the precision of the respective determi- 
nations. 



4. Discussion and Conclusions 

We have introduced a simple numerical procedure to 
evaluate the second virial coefficient of a polyatomic mole- 
cule assuming the electrostatic moments are placed at the 
center of the molecule. 

Our method should be compared, in general, to the 
expansion procedure since both techniques are based on the 
same approximations of the nonspherical intermolecular 
interaction. The expansion procedure treats the angle de- 
pendent terms of the potential as perturbations around the 
spherical contribution. Hence one obtains an expansion 
involving products of integrals of the spherical contribution 
and the electrostatic moments. For example, we have used 
eq (11) with the/j. integrals. The advantage of this approach 
is that the J x integrals need only to be calculated once for a 
given spherical potential and one has some flexibility in 
adding various angle dependent terms. 

In contrast, our method has the advantage that the / 
different angle functions and their weights [e.g., eq (7) and 
eq (10)] only have to be recorded once for any type of angle 
dependence in the potential. The method, therefore, is 
flexible as to the spherical part of the potential. 

Obviously, this flexibility has to be balanced against the 
problem that several J x integrals [eq (8)] have to be evaluated 
at a given T*. Overall, therefore, the expansion procedure is 
faster. However, we have shown [1, 2] that calculations of 
the thermophysical properties in general, and of the second 
virial coefficient, in particular, for a gas are about as 
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sensitive to a choice of the spherical part of the potential as 
they are to the various nonspherical contributions. The net 
time required, therefore, to calculate a value of Z?* for a 
given gas, is not too dependent on which of the two methods 
is used. 

When comparing the two methods one assumes that the 
expansion procedure can yield a convergent result for the 
number of terms considered. But this is known to be a 
problem for polar molecules. In fact we have pointed out 
that, for the sample calculation at T* = 0.7, r* = 1.0 
(table 3), the expansion of equation (15) has not given a 
convergent B* if j = 4. It is less appreciated that conver- 
gence should also be checked for quadrupolar molecules, as 
was mentioned with respect to carbon dioxide in reference 
[2]. [Inspection of the expansion of equation (11) indicates 
one should anticipate difficulties if, by chance, the factor 
Q* 2 /T* is close to unity.] 

As a matter of interest, we checked the convergence of 
the expansion of eq (15) for a molecule with a reduced 
quadrupole moment (squared) of 0.7 (to correspond approxi- 
mately with carbon dioxide). The reduced second virial 
coefficients from our numerical procedure were taken to be 
correct. Table 4 lists B* (numerical) at four temperatures 

TABLE 4. Comparison between values of B* for a quadrupolar molecule 

with 6* 2 = 0. 7 calculated numerically, and calculated from the expansion 

of equation (11) to various orders in 1/7 1 * 



f* 



B* 



B* [equation (11)] 



(numerical) 



0(l/7 7 * 3 ) 



0(l/r* 4 ) 



o(i/r* 5 ) 



0.7 


-6.718 


-6.2941 


-6.7729 


-6.6173 


1.0 


-3.232 


-3.1461 


-3.2434 


-3.2194 


2.0 


-.7599 


-0.7540 


-0.7605 


-0.7596 


20.0 


.5231 


.52313 


.52312 


.52312 



and the value of B* calculated from the expansion of eq (11) 
to order 1/T* 3 , 1/7 1 * 4 , and 1/7 1 * 5 , respectively. The results 
suggest that the expansion technique can fail to give a 
convergent answer at 7"* = 0.7 (and possibly at 7 1 * = 1.0) 
even to (1/7 7 * 5 ). 

In conclusion, we have demonstrated that the numerical 
technique presented here is a viable method to compute the 
second virial coefficient of polyatomic molecules whose inter- 
action is represented by the system sketched in figure 1. 
More sophisticated techniques are listed in reference [8J. 
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